UNATTAINABILITY OF A PERTURBATION BOUND FOR 
INDEFINITE LINEAR LEAST SQUARES PROBLEMS 

JOSEPH F. GRCAR* 

Abstract. Contrary to an assumption made by Bojanczyk, Higham, and Patel [3IAM J. Matrix 
Anal. AppL, 24(4):914-931, 2003|, a perturbation bound for indefinite least square problems is capa- 
ble of arbitrarily large overestimates for all perturbations of some problems. For these problems, the 
hyperbolic QR factorization algorithm is not proved to be forward stable because the error bound 
systematically overestimates the solution error of backward stable methods. 
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1. Introduction. This paper studies a perturbation bou nd for indefinite linear 
least squares problems ( Bojanczyk. Higham. and Patell . 2003 ). The problems where 
the bound is an arbitrarily large overestimate for all perturbations are characterized, 
and examples of overestimates are constructed. 

Background material is presented in section [51 The case of perturbations just to 
A is considered in section [3] The case of perturbations to both A and b is considered 
in section m 

2. Background. This section describes a canonical form for indefinite least 
squares problems and explains the perturbation bound. 

2.1. Indefinite Linear Least Squares (ILS). The ILS problem is 

X = argmin (6 — AuYj{b — Au) 

u 

where ^ G R™^", m > n, 6 e M™, x e M", J e M™''" is a diagonal matrix with only 
±1 on the diagonal, and A^JA is positive definite. The positive definite hypothesis 

reduces the problem to solving the first-order optimality con dition, A^Jjh — Ax) = 0. 

M atrix theory for indefinite inner pr oducts is develo ped bv lGohberg. Lancaster, and Rodman^ 
(1983), and for the J inner product by iHigham (2003). Questions in control theory 
that can be expressed as indefinite linear lea st squares problems are described by 
Chandrasekaran. Golub. Gu. and Savedl ( 19981 ) . They and Xu ( 2004 ) present solution 
algorithms for the ILS problem that are backward stable. 

2.2. Generalized Singular Value Decomposition. The following matrix de- 
composition provides a standard form for ILS problems. Suppose the ILS problem is 
arranged so that J has positive and negative blocks of orders and m_, 
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Note 771+ > n becaus e A*_^_A+ — A*_A- is positi ve definite. The generalized singular 
value decomposition ( Paige and Saundersi Il981 ) then takes the form 
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m+xn: C = diag(cos(6'i) j=i,...,„) 

m_ X n: S" = diag(sin(6lj) ... „ii„|^__„}) 
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where Q-\-, Q- are square, orthogonal matrices of orders m_|_, m_, respectively, X is 
invertible. and the angles satisfy, 

J > Oi > 02 > ■■■ >9n>0 and 0^ for i > to_ . 

The upper bound on the angles assures A* J A — X*{C*C — S*S)X is positive definite. 

2.3. Perturbation Bound. Suppose A + AA and b + Ab are perturbed by e in 
a "normwise relative" sense, 

IIAAIIf <6||A||;^1 . f\\AA\\F \\Abh 

(2.1) > equivalently e ma" ' 



l|A6||2<e||b||2 ^ 

where A and b are normal i zing f actors, and suppose {A + AA)* J{A + AA) is positive 
definite. Boianczvk et al. ( 20031 p. 917, eqn. 2.8) show the change to the solution 



Ax is at most 



^''"^ < IIAII^ (^\{A'JA)-'A% + UA'JA)-^ """^ 



2 ■ 



(2.2) ll^ll^ ^ ll^ll^ 

+ ||Mi||(A*JA)-iA*||.. + 0(e2), 

where r = b — Ax is the ILS residual. Bojanczyk et al. (p. 918, line -10) "believe that 
p.2p is nearly attainable." 

3. Perturbations to A. For simplicity, this section considers perturbations 
only to A. The first subsection 13.11 uses differential calculus to examine what is 
meant by an attainable perturbation bound. Subsection 13 . 21 presents the simplest ILS 
problem, for which p.2p can be an overestimate for perturbations to A. Perturbations 
to both A and b are more complicated and are deferred to section H) 

3.1. Attainable Bounds. Taylor series are the classic perturbation formulas. 
For the ILS problem, the function f{A) = {AKjA)-'^A*Jb with b fixed gives the 
dependence of the solution on the matrix, x — f{A). The perturbed value f{A + AA) 
is approximated by a power series of order p in AA under two conditions. (1) / must 
have continuous partial derivatives of order p + 1, which is always true because / 
is algebraic. (2) / must be defined on the line segment from A to A + AA, which 
amounts to (^4 + sAA)*J{A + s AA) being p ositive definite for all < s < 1. Taylor's 
formula of order p = 1 is then (jBartlel ll964L p. 245, thm. 20.16) 



(3.1) fiA + AA) - /(A) + D'fiA) AA +niAA, A) , 

J'^(A)vec(AA) 

where D^f is the first derivative of /, and TZ is the remainder. 

The first derivative is represented by a Jacobian matrix. For specificity, the 
"vec" construction is needed to order matrix entries into column vectors so partial 
derivatives with respect to matrix entries can be placed into the Jacobian matrix. 
Thus, vec(A) is the column vector of entries Aij with (j, i) in lexicographic order, 
and J7x(vec(yl)) is the Jacobian matrix of partial derivatives of / evaluated at the 
given A. Notation for Jacobian matrices has a subscript for the name of the function 
(or the dependent variable), followed by the value for the independent variable in 
parentheses. The "vec" is omitted for brevity; thus J'x{vec{A)) is written J^x{A). 
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The remainder is bounded by ||A74|||, times a coefficient independent of A A only 
when AA is restricted in some way. Let Hk be the Hessian matrix for the fc-th entry 
of f{A). The fc-th entry of TZ is ■i(vec(AA))* ■Hfcvec(AA) where Hk is evaluated at 
A + SkAA for some < < 1, which depends on A and AA. Choose (5 > so 
{A + PYJ{A + P) is positive definite when ||P||_f < 5. Let rjk be the maximum of 
||Hfc||2 over these A+P, and let r;= (ry^+jyl^... +^2)1/2^ Thus ||7^(P, A)||2 < f|lP|||. 

Subtracting f{A) from both sides of equation p.ip and applying norms gives the 
general perturbation bound, 

(3.2) \\Ax\\2<\\MA)\\2\\AA\\f + %\\AA\\1 for IIAAIIj. <,5, 

where (5, 77 only depended on A, f. This bound is remarkable in three respects. 

First, lis. ^|) is approximately attainable in the following sense. Let the induced 
norm ||i/a;(A)||2 be attained at vec(t/), which may be scaled so that ^||C/|ji;- < 
ll»/a;(^)l|2 and \\U\\p < S. For < t < 1, the perturbation tU changes the solu- 
tion by Ax{t) = f{A + tU) - /(A), which by dSH]) equals Jx{A)yec{tU) + n{tU,A). 
Triangle inequalities then imply 

(3.3) l|Aa;(t)||2 lies between \\J,{A)\\2\\tU\\F±^tU\\l. 

Thus, AA = tU can be chosen so |jAx(t)||2 is arbitrarily close to the bound (|3.2p . 

Second, the coefficient \\J^x{A)\\2 is the smallest that can appear in bounds of the 
form iS.^) . Suppose there is a comparable bound, 

||Aa;||2 < a\\AA\\F + |||AA||| for \\AA\\f <6. 

If AA = tU , then the bounds are quadratic polynomials of t. This upper bound 
exceeds the lower bound p.3|) for all < t < 1 if and only if a > || J72:(A)||2. 

Third, \\J'x{A)\\2 properly scaled gives the condition number of x with respect to 
A. Bound p.2p can be rewritten as 

(3.4) 




There may be many ways to define condition numbers, but since the scaled coefficient 
in (|3.4p is the smallest possible, the value of this coefficient must be the condition 
number of x with respect to A and the scaled norms, Xx{A). Jacobian matrices were 
first used to study conditioning by I Rice (Il966l p. 292, thm. 4). 



It is an open question to find a simple formula for the Xx{A) of general ILS 
problems. There is a formula for Jx{A) which leads to a bound for Xx(^) that 
appears in (|2.2I) . If F{A,x) — A*J{b — Ax) = (optiniality condition), then by 
implicit differentiation Jx{A) — —{Jf{x))~'^Jf{A). In the present case, by inspection 
Jf{x) = —A* J A, and by the product rule for differentiation 



Jf{A) = 



[Jrf 

{JrY 



[xiA^J ■■■ XnAKj] ^ I ®{JrY - x^ ® A^J , 



where r — b — Ax, and (g) is the Kronecker product {H (g) K has blocks consisting of K 
scaled by the corresponding entries of H). Combining these matrices by the implicit 
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differentiation formula givesQ 

(3.5) J^{A) = [A'JA)-^ ® {Jrf - x* ® {A*JA)-^AKJ . 

^ ^ ^ ^ ^ ^ 

Ml M2 
The norms of Ali and A^2 have simple formulas, which make a bound for XxiA), 

\\Mih - \\{A'JA)-^h IIHI2 IIX2II2 = Ik-lb \\{A'JA)-^A% 

(3-6) IIAIIp IIAIIp 

= IIXi - M2II2 < (IIM1II2 + IIX2II2) = 

The first line of Bojanczyk et al.'s bound (|2.2p is ;Se. 

3.2. First Example. With the choice ||b||2 = 0, ((27T]) becomes Ab = and 
e = ||AA||F/||A||i;^, and the first line of is 



(3.7) 



Il^lb \\M\f 



Comparing (|3.4p and p.7p . whether (|2.2p is nearly attainable entails the question of 
whether the formula B of p.6p is a good estimate for the condition number Xx{A). 
The answer is "no" for the problem with m = 2 and n = 1, 



02,1 



J 



The solution is 



[xi] = {A'JA)-^A*Jb 



1 



1 



ai,i + 02,1 

and the Jacobian matrix computed by differentiating x is 



Jx{A) 



dxi dxi 



dai^i' da2,i 



-1 



A 



(ai,i + ^2,1)^' (01,1+02,1)^ 



Consider matrices with the following generalized singular value decomposition, 



ai,i 




■cos(6')" 


02,1 




sin(6l) 



Choose ||A||i? = \\A\\f for the remaining scale factor, li 9 = j — a and < a < 0.3, 
then from XxiA) = sec(a) « 1.0, while from dSJ]) B = 2csc(2a) « l.O/a, in 

which 1.0 indicates values rounded to two decimal digits for these a. The formulas 



have been obtained with Mathematica (jWolfrard . 120031 ) , or they are easily calculated 
by hand. Thus, B ~ 1/a is not a good approximation to Xx{A) ~ 1. For example, if 
a = 0.01, then Xx{A) = 1.00005 and B = 100.007 rounded to six digits. This example 
has considerable cancellation in the formula p.Sp for Jx{A), 

\\Jx{A)h = ||Xi - X2II2 = 2-1/2 sec2(a) ^ 3-1/2 ^ 
WMih = II7W2II2 = 2-1/2 ^gp(2a) sec(a) « 2-3/2/a . 



^Equation 1 13. 5t was derived bv lBoianczvk et al.l | |2003|. p. 918, line 13) using perturbation algebra. 
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4. Perturbations to A and b. This section sliows tliat (|2.2I) is not nearly 
attainable in general for perturbations to both A and b. The first subsection 14. f I 
examines the types of bounds that are attainable for problems that depend on a 
matrix and a vector. Subsection 14.21 characterizes the situations in which (|2.2I) is an 
overestimate, then subsection 14.31 presents an example. 



4.1. Attainable Bounds Redux. An attainable perturbation bound for chan- 
ges only to b can be constructed as in section 13.11 The partial derivatives of a; = 
(A* JA)~^A* J& with respect to b are the entries of the coefficient matrix {A*JA)'^-^A*J, 
which is then Jx^b). Thus, like equation (|3.4p for A, so for 6, 

,4.1, 1^ < » II mi „H„e also = . 



fI|2 VI|2;||2 



Note for future reference that has a formula in terms of the matrix A^2 in the 

formulas p.6p for and B. There is no higher-order remainder term because x 

is a linear function of b. As a result, this bound can be exactly attained. 

When both A and b can be perturbed, then there are several perturbation bounds. 
The solution x — g(^A, b) = {A*^ J A)^^ A* Jb has a Taylor series analogous to (|3.ip . 

(4.2) g{A + AA,b + Ab) = g{A, b) + D^g{A, b) (AA, A6) + 5 , 

Jx{A,b)Yec{AA, Ab) 

where D^g is the first derivative of g, and S is the remainder. S is again a "second 
order" expression of AA, Ab, but for brevity a bound is not developed here. The 
"vec" construction is extended to vec{A, b) by placing the entries of b after A. This 
ordering makes J^x{A, b) into a 2 x 1 block matrix of the separate Jacobian matrices 
with respect to A and 5, Jx{A,b) = [Jx{A), Jx{b)]. 

The simplest way to convert (|4.2p into a bound applies the 2-norm to both sides 
and uses the triangle inequality on the right to treat each block of Jx{A, b) separately, 

(4.3) l|Ax||2 ^ ||^.(A)||2||A||^ IIA^II^ ^ 11^.(^^)112 ||b||2 IIA6II2 , 



2; 2 \\X\\2 WMf X 2 ||b||2 ||2;||2 



Xx(A) Xx{b) 
(4.4) < [xxiA)+xAb)]^+\\S\\2/\\xh. 

Like (|3.41 14.11) , bound ()4.3p also can be approximately attained. 

A more complicated way to convert (|4.2p into a bound requires a norm for 
vec{AA, Ab). Such a norm is |jvec(AA, A6)||inax = e of (|2.ip . The compatible norm 
for 2 X 1 block matrices has an upper bound, 

\\JxiA,b)\\= max ll^4A&)vec(AM,A.)||2 



< max 



AM,Av ||vec(AM, Av)||max 

^.(.4)112 l!A||f + 11^^(6)11, Ilb|l2^ 



AM,Av ||vec(AM, Aw) 

< |lJ.(A)|l2|lA|lp-|-|lJ.(6)||2||b||2. 
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Applying norms to (j4.2p and using this matrix norm and its bound gives two bounds^ 

\\Axh ^ \\MA,b)\\ \\sh ^ \\JAm2\\Ay + \\Mb)h\\hh \\sh 




The first of these bounds can be approximately attained in the sense of 1 

Altogether (1131 and with give four bounds with e as in (PTTj) . 




in which the 0{e^) remainder term |jS'||2/||a;||2 is omitted for clarity. The two smallest 
{i, a) can be approximately attained, but they are incommensurate. If the separate 
condition numbers Xx(A), Xx{b) are of different sizes, then the scaled perturbations 
A^, Ab can be chosen to make the sum of products (i) much smaller than the product 
of the sum and maximum (Hi). However, the bound {ii) with one condition number 
is always close to the bound (m), becaus^ 

, . XxiA)+Xxib) 

- xAAb) 

Thus, for some perturbations, bound (z) without e may be far smaller than all the 
other bounds with e. 

4.2. When ()2.2p is not Attainable. For perturbations only to A, the discus- 
sion of section [3. 1|) shows a bound is not attainable (to within high order terms) only 
if the coefficient in the bound overestimates the condition number, and p.6p shows 
^ > Xx{A) is caused by cancellation in \\M.i ~ M.2\\2 compared to ||A^i||2 + ||A^2||2- 

For perturbations to both A and b, (|4.6p shows (|2.2p is not attainable but it 
is close to an attained bound unless B + Xx{b) overestimates Xx{A) + Xx{b)- This 
situtation is much more complicated because Xx{b) my be comparable to or larger 
than the other terms. Theorem 14.31 is a precise statement of when an overestimate 
occurs. Some preparation is needed to prove the theorem. 

Lemma 4.1. If mi and "012 are positive, then 

2 .(q) mi+m2 (M 2 \mi - 1712] 
<^ ' <^ ' where ^ — 



1 + /i TO2 1 ^ fTT-i + "^2 



Proof. Taking reciprocals, multiplying by 2(mi -I-TO2), and subtracting mi +m2, 
shows the inequalities are equivalent to |mi — m2| > TO2 — nii > ~\mi — 7712]. 
Lemma 4.2. Ifu, v, w are positive and 

p= >1, then - + 1 >^"^ /9>('') -r ( - + 1 and - >'^^ p >^'^^ -- , 

V + u u 2 \ u / V 2 V 

^Boianczvk et al.l ||2003| . p. 918, line -15) use the sym b ol ■0 = Xx{A) + Xx(b)- 
^These i nequalities are remarked bv lBoianczvk et al.l l|2003l . p. 918, line -13) and are proved by 
iGrcaj 1I2OIA p. 2937, eqn. 2.10). 
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where at least one of (h) and (d) is true, but not necessarily both. 

Proof, [a) p < {w+u)/u — {w/u) + l. (c) w+u ^ pv+pu so w — pv = {p—l)u > 0, 
hence w > pv, and then w/v > p. This is the only use of p > 1. (b) li u > v, then 
p>iw + u)/{2u) = {{w/u) + l)/2. (d) liv>u, then p>{w + u)/{2v) > w/(2v). □ 

Theorem 4.3. // 

B + xAb) ^ ||.Mi||2 + ||A^2||2 ^ I|A||f||x||2 

^ xAA)+xAb) ' \\M1~M2h l|b||2 

where B, Mi, Ai2, Xx{A), Xx{b) are given by ^3.(k [.^.ip and no terms vanish, then 



2p>*Xi>p and p+i>*A2>f-l 



in which at least one of the inequalities (*) is true, but not necessarily both. Conse- 
quently, bound 112. 2\) is a uniformly large overestimate for all AA and Ab if and only 
if both Ai and A2 are much larger than 1, in which case the factor of overestimation 
is at least either Ai2~^ or A2 — 2~^. 

Proof. From (|3.6p note p> 1. With (|4.ip and some rearrangement, 

^ (IIA^lIb + ||.M2||2) + (||.M2||2 ||b||2)/(||A||f ||x||2) ^ w + u 
^ (||Xl-X2||2) + (||A^2||2||b||2)/(||A||F||x||2) V + U' 

where 



^i-(||M2||2||b||2)/(||A||f ||x||2) IIM1-X2II2 W= IIM1II2 + IIX2II2. 

Part 1, inequahty Ai > p: Note Ai = w/v, and from lemma (c) w/v > p. 

Part 2, inequality A2 > f — l^ Let mi — ||A^i||2, and let p = \mi — 7712! /(mi+m2). 
From one of the triangle inequalities, 1/Ai > p; therefore from part 1 of this proof, 
(2.i) 1/p > p. From lemma l4Tl (b). 2/(1 — p) > (mi +m2)/m2, which combines with 
(2.i) to give, (2.ii) 2/(1 — p^^) > (mi + m2)/m2. Multiplying (2.ii) by A2 and noting 
A2(mi + m2)/m2 — w/u implies, (2.iii) 2A2/(1 — p"^) > w/u. Lemma (c) is (2.iv) 
w/u > p — 1. The last two inequalities (2.iii, 2.iv) combine to 2A2/(1 — p^^) > p — 1, 
which weakens to A2 > § — 1. 

Part 3, inequality 2p >* Ai: If lemma l4T2l (d) is true, then p > {w/v)/2 where as 
noted, w/v — Ai. Rearranging gives 2p > Ai. 

Part 4, inequality p+ >* A2: If lemma (h) is true, then 2p— 1 > w/u where 
again w/u — A2(mi + m2)/m2 so altogether, (4.i) 2p — 1 > A2(mi + TO2)/m.2. From 
lemma H7T] (a) (mi + m2)/ni2 > 2/(1 + p), which combines with (2.i) to give (4.ii) 
(mi +m2)/m2 > 2/(1 + p~^). Multiplying (4.ii) by A2 and combining with (4.i) gives 
2p — 1 > 2A2/ (1 + p^^)- The latter rearranges to p + (1 — p~^)/2 > A2 which weakens 
by discarding the negative term. □ 

The formulas of theorem 14.31 resemble some given bv lBoianczvk et"aL ( 2003 . p. 



917). They introduce three quantities Ei = II7W2II2 ||b|l2/||a;||2, E2 = II7W2II2 ||A||f, 
and i?3 = ||A^i||2 ||A||i?, and they comment: "the bound (|2.2p can fail to be achieved 
for some Ab and AA only if Ei < E2 ~ E3 and there is substantial cancellation in 
the expression -{A^ JA)~'^ A* JAAx + {A*JA)~^AA*Jr for all A^. We can show in 
various special cases that these circumstances cannot arise, but we have been unable 
to establish attainability of the bound (|2.2p in general." 
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4.3. Second Example. The situation predicted by theorem 14.31 does arise for 
ILS problems with the following generalized singular value decomposition for A, 



A = 



cos(f — a) 
1 
sin(f- — a) 



1 

a- 







"1" 




'1 


1 


6 = 


1 


J = 


1 






1 




-1 



where < a < -f . The same parameter is used in both columns of A for simplicity. 
Choose the customary scale factors A = ^ and h = b. The following expressions have 
again been derived with Mathematica (|| j7a;(A)||2 is beyond easy hand calculation), 



a||f 



|b||2«3i/2 ||r||2« 2-1/2 



\\iA\JA) 

ll>^l||2 



2 « (2a)-i 
2-3/2a-i 



(A*JA)-M*||2 « (2a)-i 



\M 



2||2 



2-3/2^-1 



||X||2« 2-1/2 
||^,(A)||2« 2-131/2 

Mi-X2||2«2-i3i/2. 



In these formulas, "sa" indicates the asymptotic value as a nears 0. This example 
illustrates theorem 14.31 



X.(A) « 2-i/23i/2a~i 6«a-2 
p«6-i/2a-i Ai « 21/23-1/2^- 



xAb) « 2-1/231/2^-1 
A2 « 6-i/2a-i . 



Thus, bound (|2.2p overestimates bound (iii) by a factor of {^/6a)^^ asymptotically 
as a nears 0; therefore for all A A and A&, bound (|2.2p overestimates Aa; by at 
least the same factor. For example, if a = 0.01, then Xx{A) + Xx(b) = 244.937 and 
B + Xx{b) = 10123.1 rounded to six digits. The conditions of theorem 14.31 are. 



82.66 >* 
2p 



81.66 > 41.33 
Ai p 



41.83 >' 

p + h 



40.83 > 19.66 
A2 



2 ^ 



in which the values have been rounded to two fractional digits, and only one of the 
inequalities (*) need be true. 



5. Di scussion. The example of section H31 demonstrates that the rounding error 
analysis of Boia nczvk et al.l (|2003l ) does not prove the forward stability of the hyper- 
bolic QR factorization method to solve ILS problems. Their "overall conclusion" (p. 
929, line -12 and fol.) is that the algorithm calculates x + Ax which satisfies p.2p for 
e = 0{mnu), where u is the unit roundoff. "If we make the reasonable assumption 
that the perturbation bound (|2.2p is approximately attainable, then our rounding 
error analysis has shown that the . . . [method] is forward stable." Forward stable 
means the algo rithm sat isfies the same normwise perturbation bounds as backward 
stable methods ( Highaml . i2002.. p. 130). The complicated example of section H31 shows 
(|2.2p overestimates ||Aa:|| by an arbitrarily large factor for all perturbations A A, Ab 
of some ILS problems, in which case backward stable algorithms do satisfy a bound 
with a much smaller coefficient Xx{A) + Xx{b) B + Xx{b)- 

Nevertheless, it may be said that the hyperbolic QR factorization method is 
provisonally forward stable, in the following sense. If one or both of the quantities 
Ai, A2 of theorem 14.31 are close to 1, then p.2p approximates an attainable bound, 
and the method is forward stable for that problem. At least the size of A2 is easy to 
check: note, A2 ~ 1 is sufficient, but not necessary, for (|2.2p to be nearly attainable. 
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6. Conclusion. Of greater significance than the stability of any one algorithm 
is that the simple example of section 13.21 shows the ILS problem is among those 
for which typical normwise formulas (resembling others in numerica l linear alge bra) 
are inadequate to estimate the conditioning of the problem. .Malvshev ( 20031 ) has 



identified many problems related to least squares that have plausible error bounds 
that can be far from sharp. It is an open question how to find simple formulas for the 
condition numbers of these problems. 

Acknowledgements. I thank the editor, the reviewers, and Profs. Bojanczyk 
and Higham for corrections, remarks, and suggestions that improved this paper. 
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